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Abstract 

A new version release (2.0) of the molecular simulation tool ms2 [S. Deublein et al., Comput. Phys. Commun. 
182 (2011) 2350] is presented. Version 2.0 of ms2 features a hybrid parallelization based on MPI and OpenMP 
for molecular dynamics simulation to achieve higher scalability. Furthermore, the formalism by Lustig [R. 
Lustig, Mol. Phys. 110 (2012) 3041] is implemented, allowing for a systematic sampling of Massieu potential 
derivatives in a single simulation run. Moreover, the Green-Kubo formalism is extended for the sampling of the 
electric conductivity and the residence time. To remove the restriction of the preceding version to electro-neutral 
molecules, Ewald summation is implemented to consider ionic long range interactions. Finally, the sampling of 
the radial distribution function is added. 


1. Introduction 


Molecular modeling and simulation is a technology central to many areas of research in academia and in¬ 
dustry. With the advance of computing power, the scope of application scenarios for molecular simulation is 
widening, both in terms of complexity of a given simulation and in terms of high throughput. Nowadays, e.g. 
the predictive simulation of entire phase equilibrium diagrams has become feasible. However, in order to rely on 
simulation results, the methodology needs to be sound and the implementation must be thoroughly verified. In 
its first release ll, we have introduced the molecular simulation tool ms2. Results from ms2 have been verified 
and the implementation was found to be robust and efficient. 

As described in Section^ in Version 2.0 of the simulation tool ms2 the existing molecular dynamics (MD) 
MPI parallelization was hybridized with OpenMP, leading to an improved performance on multi-core processors. 
Furthermore, the new release offers a wider scope of accessible properties. In particular, ms 2 was extended to 
calculate Massieu potential derivatives in a systematic manner, cf. section[3] This augments the range of sampled 
properties significantly and, as was demonstrated in [23, it allows to straightforwardly develop competitive 
fundamental equations of state from a combination of experimental VLE data and molecular simulation results. 
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Lastly, besides being now capable of simulating ionic substances, the time and memory demand for calculating 
transport properties was reduced significantly (section^. 

ms2 is freely available as an open source code for academic users at www. ms-2. de 

2. Hybrid MPI & OpenMP Parallelization 

The molecular simulation tool ms2 focuses on thermodynamic properties of homogeneous fluids. Therefore, 
systems investigated with ms2 typically contain on the order of 10 3 molecules. While for Monte Carlo simula¬ 
tions a perfect scaling behavior up to large numbers of cores can be trivially achieved, MD domain decomposition 
- the de facto standard for highly scalable MD - is not feasible for such system sizes, because the cut-off radius 
is in the same range as half the edge length of the simulation volume. This excludes domain decomposition and 
limits the scalability of the MPI parallelization. The present release of ms2 features an OpenMP parallelization, 
which was hybridized with MPI. At the point where MPI communication becomes a bottleneck, a single process 
still has enough load to distribute to multiple threads, improving scalability. 

Three parts of ms2 were parallelized with OpenMP: the interaction partner search, the energy and the force 
calculations. All OpenMP parallel regions rely on loop parallelism, as the compute intensive parts of the al¬ 
gorithm all feature a loop over the molecules. In the force calculation, race conditions need to be considered, 
because every calculated force is written to both interacting molecules. Introducing atomic updates or critical 
sections leads to massive overheads. Instead, it is more efficient to assign forces from individual interactions to 
the elements of a list (or an array) which is subsequently summed up. The same holds true for torques. 

In Figure[I]the speed-up of hybrid MPI/OpenMP vs. pure MPI is plotted for 2’048 cores, varying the number 
of threads per MPI process and the number of molecules in the simulation volume. As can be seen, using 2 
to 4 threads per MPI process delivers a speed-up of around 20% for 2’048 cores. The evaluation of the hybrid 
parallelization algorithm was performed on a CRAY XE6 Supercomputer at the High Performance Computing 
Center in Stuttgart, which has an overall peak performance of one PFLOPS. It consists of 3552 nodes, each 
equipped with two AMD Opteron 6276 (Interlagos) processors. Each processor has 16 cores, sharing eight FPUs 
(Floating Point Units). Nodes are equipped with 32 GB RAM and are interconnected by a high-speed CRAY 
Gemini network. Additional runtime performance comparisons with the simulation tool GROMACS | 3] are listed 
in Table Q] 


3. Massieu potential derivatives 


ms2 version 2.0 features evaluating free energy derivatives in a systematic manner, thus greatly extending 
the thermodynamic property types that can be sampled in single simulation runs. The approach is based on the 
fact that the fundamental equation of state contains the complete thermodynamic information about a system, 
which can be expressed in terms of various thermodynamic potentials |4], e.g. internal energy E(N, V, S), en¬ 
thalpy H(N,p, ,5'), Helmholtz free energy F(N. V, T) or Gibbs free energy G(N. p. T ), with number of particles 
TV, volume V , pressure p, temperature T and entropy S. These representations are equivalent in the sense that 
any other thermodynamic property is essentially a combination of derivatives of the chosen form with respect 
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Table 1: Runtime performance results with ms2 release 2.0 and GROMACS v4.6.5 [j^] for MD simulations with 
pure water at 298.15 K and 55.345 mol • dm -3 . The number of time steps were 100 000 for every simulation, the 
cutoff radius was identical for simulations with the same number of particles. All simulations were performed 
on the same computer cluster. 


cores 

threads 

N j 

gromacs / s 

ms2(RF) / s 

ms2(EW) / s 

8 

8 MPI 

500 

164 

416 

785 

8 

8 MPI 

1000 

299 

874 

1607 

8 

8 MPI 

2000 

1284 

4461 

6777 

16 

16 MPI 

500 

95 

233 

415 

16 

16 MPI 

1000 

166 

477 

848 

16 

16 MPI 

2000 

678 

2298 

3506 

32 

32 MPI 

500 

62 

152 

245 

32 

32 MPI 

1000 

106 

296 

487 

32 

32 MPI 

2000 

361 

1286 

1898 

64 

64 MPI 

500 

40 

119 

166 

64 

64 MPI 

1000 

65 

228 

324 

64 

64 MPI 

2000 

220 

814 

1261 

128 

128 MPI 

500 

38 

105 

131 

128 

128 MPI 

1000 

51 

197 

247 

128 

128 MPI 

2000 

147 

557 

727 

8 

1 MPI, 8 OMP/MPI 

500 

167 

483 


8 

1 MPI, 8 OMP / MPI 

1000 

323 

975 


8 

1 MPI, 8 OMP / MPI 

2000 

1416 

4831 


16 

2 MPI, 8 OMP / MPI 

500 

105 

253 


16 

2 MPI, 8 OMP / MPI 

1000 

186 

517 


16 

2 MPI, 8 OMP / MPI 

2000 

763 

2514 


32 

4 MPI, 8 OMP / MPI 

500 

75 

167 


32 

4 MPI, 8 OMP / MPI 

1000 

121 

316 


32 

4 MPI, 8 OMP / MPI 

2000 

418 

1362 


64 

8 MPI, 8 OMP / MPI 

500 

60 

119 


64 

8 MPI, 8 OMP / MPI 

1000 

92 

217 


64 

8 MPI, 8 OMP / MPI 

2000 

261 

785 


128 

16 MPI, 8 OMP / MPI 

500 

49 

101 


128 

16 MPI, 8 OMP / MPI 

1000 

74 

172 


128 

16 MPI, 8 OMP / MPI 

2000 

170 

496 



(N) Number of water molecules 

(RF) simulations with reaction field correction. 

(EW) simulations with Ewald summation. 
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#MPI Threads x #OpenMP Threads 

Figure 1: Speed-up of hybrid MPI/OpenMP vs. pure MPI for 2048 cores, varying number of threads per MPI 
process and 8192 molecules (solid circles), 4096 molecules (empty circles), 2048 molecules (solid triangles), 
1024 molecules (empty triangles) 


to its independent variables. The form F/T(N, V. 1/T), known as the Massieu potential, is preferred in molec¬ 


ular simulations due to practical reasons 


ilfl- 


The statistical mechanical formalism of Lustig allows for the 


simultaneous sampling of any A r mn in a single NVT ensemble simulation for a given state point 

d m+n (F/(RT)) om n _ A 

Qfjm dp n ^ P ~ J ^ mn ~ A-mn + A mn , (1) 

where R is the gas constant, (3 = 1/T and p = N/V. A mn can be separated into an ideal part A l mn and a residual 


part A r mn |9|]. The calculation of the residual part is the target of molecular simulation and the derivatives /l ' l 0 , 
.Apt, A 20 , Ah, Aq 2 , Ag 0 , A'h and A r 12 were implemented in ms2 for NVT ensemble simulations. The ideal 
part can be obtained by independent methods, e.g. from spectroscopic data or ab initio calculations. However, it 
can be shown that for any A mn = A' rnn + AL n , where n > 0, the ideal part is either zero or depends exclusively 
on the density, thus it is known by default |g|. Note that the calculation of /l ( ' l0 still requires additional concepts 
such as thermodynamic integration or particle insertion methods. From the first five derivatives Am, ,4 ( ri, .490. 
An, Aq 2 every measurable thermodynamic property can be expressed (see the supplementary material for a list 
of properties) with the exception of phase equilibria. A detailed description of the implementation is in the sup¬ 
plementary material, here, only an overview is provided. 

The calculation of the derivatives up to the order of n = 2 requires the explicit mathematical expression of 
dU/dV and d 2 U/dV 2 with respect to the applied molecular interaction pair potential and has to be determined 
analytically beforehand [Q, f)]- The general formula for d n U/dV n can be found in Ref. |8]. For common 
molecular interaction pair potentials, like the Lennard-Jones potential ifiol 11], describing repulsive and disper¬ 
sive interactions, or Coulomb’s law, describing electrostatic interactions between point charges, the analytical 
formulas for dU/dV and d 2 U/dV 2 can be obtained straightforwardly. 

As molecular simulation is currently limited to operate with considerably fewer particles than real systems, the 
effect of the small system size thus has to be counter-balanced with a contribution to U and d n U/dV n called 

nn 

long range correction (LRC) IHOU l 11- The mathematical form of the LRC depends on the molecular interac- 
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tion potential and the cut-off method (site-site or center-of-mass cut-off mode) applied. For the Lennard-Jones 


potential, the LRC scheme was well described in the literature for both the site-site [5, 


mass cut-off mode ||8j. 


1 3fl • The reaction field method 0 was the default choice in the preceding version of 


12] and the center of 


ms2 for the LRC of electrostatic interactions modelled by considering charge distributions on molecules. The 
usual implementation of the reaction field method combines the explicit and the LRC part in a single pair poten¬ 


tial fl4, 15] from which ()'" U/dV n (including the LRC contribution) is directly obtainable. However, practical 
applications show that the electrostatic LRC of dU/dV and d 2 U /dV 2 can be neglected in case of systems for 
which the reaction field method is an appropriate choice. E.g., the contribution of the electrostatic LRC for a 
liquid system ( T = 298 K and p = 45.86 mol/1) consisting of only 200 water and 50 methanol molecules with 
a very short cut-off radius of 20% of the edge length of the simulation volume is still << 1% for both DU/dV 
and d 2 U/dV 2 . The supplementary material contains detailed elaborations on the LRC for the Lennard-Jones 
potential. 


4. Algorithmic Developments 

Transport property calculations. In ms2, transport properties are determined via equilibrium MD simulations 
by means of the Green-Kubo formalism ifltil . This formalism offers a direct relationship between transport coef¬ 
ficients and the time integral of the autocorrelation function of the corresponding fluxes. An extended time step 
was defined for the calculation of the fluxes, the autocorrelation functions and their integrals. The extended time 
step is n times longer than the specified MD time step, where n is a user defined variable. The autocorrelation 
functions are hence evaluated in every n-th MD time step. As a consequence, the memory demand for the au¬ 
tocorrelation functions was reduced and the restart files, which contain the current state of the autocorrelation 
functions and time integrals, become accordingly smaller. In addition, the overall computing time of the MD 
simulation was reduced significantly. 


Ewald summation. Ewald summation 110,11] was implemented for the calculation of electrostatic interactions 
between point charges. It extends the applicability of ms2 to thermodynamic properties of e.g. ions in solutions. 
In Ewald summation, the electrostatic interactions according to Coulomb’s law are divided into two contributions: 
short-range and long-range. The short-range term includes all charge-charge interactions at distances smaller 
than the cut-off radius. The remaining contribution is calculated in Fourier space and only the final value is 
transformed back into real space. This allows for an efficient calculation of the long-range interactions between 
the charges. The algorithm is well described in literature. Currently, some of the new features, the calculation 
of Massieu potential derivatives and Hybrid MPI & OpenMP Parallelization for MD, are not available together 
with Ewald summation. 

5. Property Calculations 

Radial distribution function. The radial distribution function (RDF) g(r) is a measure for the microscopic struc¬ 
ture of matter. It is defined by the local number density around a given position within a molecule p L (r ) in 
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relation to the overall number density p = N/V 

9 ^ p p dV 4tt r 2 p dr ' ( } 

Therein, dJV(r) is the differential number of molecules in a spherical shell volume element dV, which has the 
width dr and is located at the distance r from the regarded position. g(r) can be evaluated for every molecule of 
a given species. 

In the present release of ms2, the RDF can be calculated during MD simulation runs for pure components and 
mixtures on the fly. The RDF is sampled between all LJ sites. In order to evaluate RDFs for arbitrary positions, 
say point charge sites, superimposed dummy LJ sites with the parameters a = e = 0 have to be introduced in the 
potential model file by the user. 


Electric conductivity. The evaluation of the electric conductivity a was implemented in ms2 version 2.0, being 


a measure for the flow of ions in solution. The Green-Kubo formalism 
(T and the time-autocorrelation function of the electric current flux j,At) 

1 


16] offers a direct relationship between 


<7 = 


(je(t) ■ j e (0))dt , 


(3) 


3Vk B T J 0 

where fc B is Boltzmanns constant. The electric current flux is defined by the charge (p : of ion k and its velocity 
vector Vk according to 

je(t) = ^2qk-v k {t), (4) 


fc=l 


where Nj is the number of molecules of component j in solution. Note that only the ions in the solution have to 
be considered, not the electro-neutral molecules. For better statistics, a is sampled over all independent spatial 
elements of j e (t). 


Thermal conductivity of mixtures. In the previous version of ms2 the determination of the thermal conductivity 
by means of the Green-Kubo formalism was implemented for pure substances only. In the present release, the 
calculation of the thermal conductivity was extended to multi-component mixtures. The thermal conductivity A 
is given by the autocorrelation function of the elements of the microscopic heat flow J* 


A = 


Vk B T 2 


dt(r q (f)-j*{ 0 )). 


(5) 


In mixtures, energy transport and diffusion occur in a coupled manner, thus, the heat flow for a mixture of n 
components is given by 111811 


n Ni 


1 ^ 
j, = Jee 


2=1 k —1 




m\ (v?) + wflfwf + E u ( r fJ) 

2=1 


n n Ni Nj 

tEEEE 7 E( v ' • 


Ni 


2=1 j =1 k— 1 If^k 


f-^nnkl 

u i j 


■w?rg)-5>E v ?- 


(6) 


2=1 k =1 
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where wf is the angular velocity vector of molecule k of component i and if its matrix of angular momentum 
of inertia, u (Yff) is the intermolecular potential energy and rfj is the torque due to the interaction of molecules 
k and /. The indices i and j denote the components of the mixture, hi is the partial molar enthalpy. It has to be 
specified as an input in the ms2 parameter file and can be calculated from NpT simulations. 


Residence time. The residence time r 7 - defines the average time span that a molecule of component j remains 
within a given distance r,; 7 around a specific molecule i. It is given by the autocorrelation function 


1 


j(0) 


t=o \ n «i(0) 


^ 0fe(f)Bfc(O) ) dt , 


(7) 


k =1 


where t is the time, n,; 7 - (0) the solvation number around molecule i at t = 0 and 0 is the Heaviside function, 
which yields unity, if the two molecules are within the given distance, and zero otherwise. Following the proposal 
of Impey et al. [19|], the residence time explicitly allows for short time periods during which the distance between 
the two molecules exceeds r, 7 . Also, the solvation number n, :) can be evaluated on the fly 


rtij 



Jo 

where pj is the number density of component j and r tnill 
calculated. 


r 2 gij(r)dr, (8) 

is the distance up to which the solvation number is 
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